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Exploring supersolidity in naturally occurring and artificially designed systems has been and will 
continue to be an area of immense interest. Here, we study how superfluid and charge-density- 
wave (CDW) states cooperate or compete in a minimal model for hard-core-bosons (HCBs) coupled 
locally to optical phonons: a two-dimensional Bose-Holstein model. Our study is restricted to the 
parameter regimes of strong HCB-phonon coupling and non-adiabaticity. We use Quantum Monte 
Carlo simulation (involving stochastic-series-expansion technique) to study phase transitions and to 
investigate whether we have homogeneous or phase-separated coexistence. The effective Hamiltonian 
involves, besides a nearest-neighbor hopping and a nearest-neighbor repulsion, sizeable double¬ 
hopping terms (obtained from second-order perturbation). At densities not far from half-filling, 
in the parameter regime where the double-hopping terms are non-negligible (negligible) compared 
to the nearest-neighbor hopping, we get checkerboard-supersolidity (phase separation) with CDW 
being characterized by ordering wavevector Q = (rr^). 

PACS numbers: 71.38.-k, 67.80.kb, 71.45.Lr, 74.20.Mn 


I. INTRODUCTION 

Whether diagonal long range orders [such as charge- 
density-wave (CDW) and spin-density-wave (SDW)] and 
off-diagonal long range orders [such as superconduct¬ 
ing and superfluid (SF) states] can coexist homoge¬ 
neously in correlated electronic systems is a central is¬ 
sue in condensed matter physics. Coexistence of su¬ 
perconductivity and CDW has been studied in three- 
dimensional systems^ (such as BaBiCU doped with 
K or Pb), quasi-two-dimensional systems^ (such as 
the dichalcogenide 2H — TaSe 2 and NbSe 2 ) and quasi- 
one-dimensional systems^ (such as the trichalcogenide 
NbSe 3 and doped spin ladder Sri 4 Cu 24 C> 4 i). 

Supersolidity, defined as homogeneous coexistence of 
superfluidity and crystalline order, was theoretically pro¬ 
posed more than 4 decades ago^ and has been a subject of 
debate since then. The first experimental clairn^ of ob¬ 
serving supersolidity in helium-4 further intensified the 
debate and enhanced the interest in understanding the 
phenomena. The general consensus^ is that a supersolid 
(SS) is not realizable in a perfect hep crystal of 4 He; nev¬ 
ertheless, superflow can occur along vacancies that can 
collect near extended structural defects^. Thus, the oc¬ 
currence of supersolidity in bulk solid helium-4 is being 
seriously questioned. 

Cold-atom systems offer another opportunity for re¬ 
alization of supersolidity. Theoreticall y 3 : 4 : 12 — , lattice 
bosons with various types of interactions in diverse ge¬ 
ometries have yielded supersolidity. However, there has 
been no experimental creation of optical lattices with ef¬ 
fective long-range interactions that produce supersolid¬ 
ity. Furthermore, experimental techniques to detect sig¬ 
natures of supersolidity also need to be developed for 
optical lattices^. 

There have been numerous studies of supersolidity 
involving hard-core-bosons (HCBs) ^ 12 ’ 13 . A lattice 
model for quantum liquids, such as the interacting 


bosonic helium-4 at low temperatures, needs a hard-core 
constraint to account for the exclusion of occupation of 
more than one atom at each lattice poin t 18 : 19 . The be¬ 
havior of the ground state and low temperature excita¬ 
tions of such systems are largely controlled by couplings 
with phonons. Furthermore, local Cooper pairs [compris¬ 
ing of two electrons (of opposite spin) at a site] can also 
be regarded as HCBs^S. In Bismuthates, such HCBs cou¬ 
ple to the breathing mode of the oxygen cage surrounding 
the Bismuth ion s 21 : 22 . 

Here, in this article, we study a two-dimensional (2D) 
Bose-Holstein (BH) model for HCBs on a square lattice 
where they can hop to nearest-neighbor (NN) sites and 
experience the HCB-phonon interactions via a Holstein- 
type term. Previously, exact diagonalization calcula¬ 
tions were done on this model^ for a small system (i.e., 
4x4 lattice) to study the resulting phase diagram. Here 
we use stochastic-series-expansion (SSE) based quantum 
Monte Carlo (QMC) technique to simulate large-size lat¬ 
tices so that various phases in the thermodynamic limit 
can be identified more clearly. Unlike the t — V model, 
a SS is realized in our BH system due to non-negligible 
transport within the same sublattice. At densities not 
far from half-filling and at sufficiently large HCB-phonon 
couplings, phase coexistence occurs; furthermore, in the 
phase-coexistence region, the system tends to phase sep¬ 
arate at stronger couplings. 

Our paper is organized as follows: section II deals with 
a discussion of the BH Hamiltonian, its transformations 
and its mapping to the equivalent spin model. Section 
III covers a description of the numerical method and the 
observables employed to characterize the orderings. In 
section IV, we detail our results and the corresponding 
analysis, both with and without the presence of same- 
sublattice hopping terms. Lastly, in section V, we sum¬ 
marize our results and draw conclusions. 
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II. FORMULATION 


The BH Hamiltonian is given by 

H = -t^b^bj+s + w 0 E a ] a i + guj 0 Y n ^ a 3 + aj),(l) 
j,s i i 

where cq and bi denote the annihilation operators for 
phonons and HCB particles, respectively, and rq (= b\bi) 
is the number operator for HCBs at site i. Furthermore, 
t is the amplitude for hopping at NN sites denoted by 
6 and ojq is the frequency of optical phonons^i. An ef¬ 
fective Hamiltonian for the HCB particles is obtained 
by first transforming this BH Hamiltonian to the pola- 
ronic frame of reference (using the Lang-Firsov trans¬ 
formation) and then performing perturbation theory as 
detailed in Refs. [23] and [25|. Interestingly, second-order 
perturbation theory yields a two-step hopping which pro¬ 
duces, in a 2D square lattice, both next-nearest-neighbor 
(NNN) and next-to-next-nearest neighbor (NNNN) hop¬ 
ping terms besides the usual NN hopping term in the 
Hamiltonian. Moreover, a NN repulsion also results from 
two-step virtual hopping back and forth. 

So we get an effective t\—t 2 — t^ — V Hamiltonian for 
HCB particles on a 2D square lattice^: 


He 


-g 2 u o E n i ~ ^ E b j b i+s -faE b jbj+8' 

3 3,8 3,8' 

1 3 E h ) b 3+8" n j+«s)> ( 2 ) 

3,8" 3,8 


where S' and 5" denote NNN and NNNN sites respec¬ 
tively; fi = texp(-g 2 ), t 2 = (2ff/w 0 )/i(ff), £3 = £ 2/2 and 
V = (tf/cj 0 )[4f 1 (g)+2f 2 (g)]; /1 (g) = 9 2n /{n\n) and 

f 2 (g) = X^m=i g 2< ' n+m ' > /[n\m\(n + m)\. In the regime 

g 2 

g > 1 , we can make the approximations fi(g) ~ ^ 3 - and 
2fl 2 

[/2 (< 7 ) + 2/i (g)] ~ with the approximations becom¬ 
ing exact for g —> 00. The small parameter is t/(guJo) and 
is obtained from [U/wo] 1 ^ 2 (see Ref. [ 2 ^ for details); our 
perturbation analysis is done in the nonadiabatic regime 
( t < u) 0 ) and at strong coupling ( g > 1). 

Since all the hoppings are non-frustrated, a QMC simu¬ 
lation of the system does not suffer from the negative sign 
problem. We employ SSE technique in our QMC simu¬ 
lation and investigate the co-existence or competition of 
CDW and superfluidity in various regimes of the parame¬ 
ter space. Here we should mention that our t\—t 2 —to — V 
model for HCBs is equivalent to an extended XXZ spin- 
1/2 Hamiltonian as shown below: 


H=Y [Ji*S?S? + %( S+Sr + H.c.)]+ 


<ij> 

J'lxy 
2 

J 3 xy 


E (StST + U. c.)+ 


«i,j» 


E (StST+R.c.)-hoY S ^ ( 3 ) 


«<l,3' 


TABLE I. Ji z and J 2xy in terms of J\ xy [in Eq. l[3|)] at various 
values of the HCB-phonon coupling g and at t/uj 0 = 1. 


g 

0.5 

1.0 

1.5 

2.0 

2.5 

3.0 

J lz 

0.444 

1.355 

2.725 

8.017 

45.485 

478.571 

J2xy 

0.415 

0.970 

0.960 

0.647 

0.395 

0.255 


where < i,j >, << i,j >>, and <<< i,j >>> stand 
for NN, NNN, and NNNN pairs respectively. Further¬ 
more, the operators for the HCBs are related to those 
of spin-1/2 particles as: 5| = rij — \ and S+ = by A 
comparison of the parameters in Eqs. © and m yields: 
Jiz = V , J\ X y — 2ti, J 2 xy = 2Q, J'ixy — 2^3 and 

ho = g 2 u}g. Now, the magnetization of the system can be 
tuned by using an external magnetic field; then, a term 
—hJi X y JT Sf should be added to the Hamiltonian in Eq. 
© where the magnetic field h is given in units of J\ xy . 

Presence of hopping terms for HCBs indicates that su¬ 
perfluidity (i.e., spontaneous breaking of the global U(l) 
gauge symmetry) can exist in the system. On the other 
hand, a large interaction strength suggests the possibility 
of a CDW. Thus, our objective is to study the compat¬ 
ibility of these two long range orders. Now, these two 
orders can coexist either in a phase separated form or 
homogeneously as a SS. It should be pointed out that, a 
t — V model on a square lattice does not show a thermo¬ 
dynamically stable SS phase for HCBs^. On the other 
hand, striped SS behavior is found away from half filling 
when NNN repulsion (V 2 ) is considered in the t — \\—V 2 
model^— —. 


III. NUMERICAL CALCULATIONS USING 
SSE-QMC 

We now give details of the SSE-QMC simulation of 
our t\ — t 2 — to — V model, or equivalently, our ex¬ 
tended XXZ spin-1/2 model. Finding the phase dia¬ 
gram in the present problem requires exploring various 
limits of the parameters (including high anisotropy in 
our spin model). In our numerical computations, we 
used directed loop update for efficient sampling of the 
configuration s 28 ! 29 . The ground state properties are cap¬ 
tured by simulating at low enough temperatures, i.e., 
/? ~ L with L being the linear dimension of the square 
latticed. We employ /3 = 3L/2 since our calculations 
with (3 = 2L yield the same values for the observables 
(within the error bars of the calculations). From calcu¬ 
lations involving various large system sizes, we infer the 
results in the thermodynamic limit. 

As can be seen from the expressions of the two-spin ma¬ 
trix elements for Heisenberg spin models used in various 
SSE-QMC studie s 28 ! 29 , a positive parameter e is intro¬ 
duced to ensure the positivity of all the matrix elements 
(see Appendix A for details). This is necessary so that 
they can be treated as probabilities. The value of e can 
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also affect the autocorrelation time of the desired vari¬ 
able. We found that keeping the numerical value of e 
equal to at least a quarter of the anisotropy parameter 
Jiz/Jixy, particularly near the transition region, helps 
keep the data in various bins uncorrelated when each bin 
contains data from large number of Monte Carlo sweeps 
(i.e., at least 15,00,000); for the value of e used in the 
XXZ model, see Ref. [28|. 

In this work we are concerned with the diagonal order 
parameter [i.e, the structure factor at the Neel 

ordering vector Q = (7r, 7 t)] and the off-diagonal order 
parameter of the SF density ps■ A general expression for 
S(Q) (for our HCB system) is given as 

S(Q) = ~ ( 4 ) 

where () denotes the ensemble average. Since rii (or 
Sf in our spin model) are diagonal in the basis, a QMC 
average can be computed easily. 

The SF density^! is given by ps = 1 /N {d 2 F / dO 2 ') g _ Q , 
where F is the free energy in the presence of twisted 
boundary conditions with angle of twist 0. This is an 
off-diagonal order parameter. In a QMC calculation, the 
SF density along x-direction is calculated using ps x =< 
(iV+ — N~) 2 > //3N where ./V+ and N~ represent the 
total no. of Hamiltonian operators transporting spin in 
the positive and negative x-directions, respectively^!. 

A few benchmarking comparisons of our SSE results 
with exact results for a J\ — J 2 model are shown in the 
supplementary material. 

IV. ANALYSIS OF RESULTS AND 
TECHNICALITIES 

Our calculations for the 2D model aim at 

obtaining the ground state phase diagram of the system 
and also at extending (to the thermodynamic limit) the 
findings presented by a Lanczos study on a small cluster 
in Ref. |23|. We consider t/uo = 1 in the present study. 

To obtain the phase diagram for our system, the in¬ 
terplay between the diagonal and off-diagonal orders and 
the transition between them needs thorough investiga¬ 
tion. As antiferromagnetic order breaks SU(2) symme¬ 
try whereas a SF phase breaks U(l) symmetry, a phase 
transition between them (according to Landau theory) 
cannot be of second-order type. Furthermore, a transi¬ 
tion to a phase separated state occurs when the system 
undergoes a first-order transition. 

We work in the grand canonical ensemble where there 
is no constraint on the HCB particle number (or the mag¬ 
netization in the equivalent extended XXZ model). In 
Fig. [T] we show the variation of m z , S(n, 7 r), and ps with 
magnetic field h (expressed in units of J\ xy ) for g=1.75 
for both our extended model and its XXZ version (i.e., 
J 2 xy = J%xy = 0) in L x L square lattices with L = 10, 14 
and 16. The comparison shows that the results for differ¬ 
ent size lattices almost coincide. At small values of the 



FIG. 1. (Color online) Identifying sufficiently large sys¬ 
tem sizes to capture quantum phase transition, ps, 

S (71", 7 r ) and m z vs magnetic field h in L x L square lat¬ 
tices with L = 10 (line), 14 (circle) and 16 (triangle up) for 
t/uio = 1 and g = 1.75 when (a) J 2 xy = ‘ZJzxy 7 ^ 0 and (b) 

Jlxy — J^xy — 0 . 


magnetic field h, the system manifests half-filling in the 
case of HCBs (or zero magnetization in the case of the 
equivalent spin model); the CDW phase is formed with 
maximum values of S(tt,tt) with SF density simultane¬ 
ously assuming zero value. At large values of h, before 
the system is at complete filling, S(n, 7r) decreases to zero 
while ps becomes finite manifesting the SF phase. In 
the intermediate magnetic field region, phase transition 
occurs with both the orders coexisting in our extended 
XXZ model. The plot of magnetization m z [or (p — 1/2) 
where p is the density] in Fig. QJb), indicates a discrete 
jump during the transition in the case of the XXZ model. 
This is a typical signature for a first-order transition. 
Contrastingly, a continuous smooth increase in m z is ob¬ 
served for the extended XXZ model clearly ruling out 
the possibility of a phase separation. Calculations us¬ 
ing a canonical ensemble shows an inhomogeneous phase 
coexistence due to phase separation (PS)^. 

In the region adjacent to m z = 0 in Fig. [l](a), where 
overlap of and ps are observed, the system dis¬ 

plays a SS phase. Analysis of various large system sizes 
confirms the picture in Fig. [2a) that there is indeed 
a homogeneous coexistence of CDW and SF long-range 
orders in the system. 

Next, we will do a further detailed case-by-case study 
for the system with and without the effects of same- 
sublattice (i.e., NNN and NNNN) hoppings. 


A. Considering only NN hopping (t 2 = £3 = 0) 

Here we study the case of £2 = £3 = 0 in Eq. m, 
i.e., the bare XXZ model. Calculations are done on a 
16 x 16 lattice. We find that the system loses its CDW 
order at half-filling for values of the coupling g below the 
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m z (p - 1/2) 


FIG. 2. (Color online) Quantum phase diagrams 
at various magnetizations (fillings) and adiabaticity 
t/cjo = 1.0. The calculations are for a 16 x 16 lattice and for 
our BH system [using Eq. ([2)>] by (a) considering t 2 = t 3 = 0 
and (b) including all the interactions and hoppings. 


critical value of g c = 0.82 corresponding to the Heisen¬ 
berg point of the XXZ model. For smaller values of g, 
superfluidity develops for all values of filling between 0 
and 1. The phase diagram for the XXZ model is depicted 
in Fig. ©a). Our calculated phase diagram of the XXZ 
model is compatible with Fig. 1, Fig. 2a and Fig. 2b 
of Ref. [Hi. At half filling (i.e., p = 1/2 or m z = 0), the 
Heisenberg point denotes the boundary between CDW 
and SF phases. 

Fig. Oja) shows the jumps in magnetization m z as soon 
as we increase g beyond g c . The magnitude of the jump 
increases as the coupling g increases. This jump implies 
a first-order transition and indicates a phase separated 
coexistence of the CDW and SF phases. A histogram 
analysis can also capture the jump in m z values. In Fig. 
©b), plotted for g = 1.50, the m z histograms change 
when magnetic field values are varied. For instance, when 
the magnetic field is set at h = 4.44, m z = 0 even when 
large number of Monte Carlo sweeps were used in our 
QMC computation. On increasing the magnetic field to 
h = 4.46, the m z values start showing a double-peaked 
structure (with one peak at m z = 0 and another peak 
at m z « 0.18) which is indicative of phase separation 
[see also inset in Fig. ©b)[. Then, at a slightly higher 
magnetic field of h = 4.48, all the m z values seem to 


be centered around a mean value close to 0.18. This 
manifests the first-order transition. All the discontinuous 
transitions are due to inhomogeneous coexistence of the 
CDW and SF phases. 


B. Considering all hoppings 

The phase diagram for our BH model [obtained from 
Eq. (©)] is depicted in Fig. ©b). On including the effects 
of NNN and NNNN hoppings (i.e., = 2 f 3 ^ 0) in 

the t\ — t 2 — £3 — V model of Eq. there can be a 

difference in the densities in the two sublattices owing to 
NN repulsion and same-sublattice hopping. As shown in 
Table. E at intermediate values of g (i.e., g ~ 1 ), NNN 
hopping is comparable to NN hopping; consequently, a 
SS state can occur. Whereas at larger values of g (i.e., 
for g > 2.5), NNN hopping is fairly smaller than NN 
hopping and we can expect the same behavior as in the 
t — V (or the XXZ) model. On account of particle-hole 
symmetry in our model, the phase diagram is symmetric 
about half-filling. 

We will now examine various features of the phase di¬ 
agram of Fig. ©b). The phase diagram for our model 
was obtained by identifying the transition regions and 
the nature of the various phases. The variation of mag¬ 
netization with magnetic field in a 16 x 16 lattice is 
shown in Fig. [ 6 ] in the transition region. The results 
for g = 1.5, 1.75, & 2.0 show that the magnetization 
increases gradually without any jump as the magnetic 
field is increased. Hence, in Fig. ©a), for g = 1.75 & 2.0 
CDW and SF phases coexist homogeneously resulting in 
a SS state. Here we must mention that Fig. ©a) (de¬ 
picting simultaneous coexistence of CDW and SF phases 
through non-zero values of and ps) corroborates 

this conclusion. 



FIG. 3. (Color online) Magnetization plots showing dis¬ 
continuous transitions in the XXZ model, (a) Magneti¬ 
zation m z vs magnetic field h in a 16 x 16 lattice at /3 = 1.575, 
at t/u>o = 1, and for different values of g. (b) Magnetiza¬ 
tion histograms in the transition region at g = 1.5; two peak 
structure of the histogram for h = 4.46 is highlighted in the 
inset. 
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FIG. 4. (Color online) Determining cut-off for the 
anisotropy parameter A = Jiz/Jixy when J2xy 7^ 0. 

Results of (a) m z vs h and (b) S(n,n) vs m z for a 16 x 16 
lattice, for t/uio = 1, and for g = 2.5 but with A = 10, 15 
and 20. 


As a matter of fact, large anisotropy [i.e., large values 
of V/ti in Eq. (O] in the model requires large simulation 
time. Thus, as we increase the value of g and thereby 
increase the value of V or J\ z (see Table HJ, the numeri¬ 
cal calculation suffers from appreciable slowing down and 
with our computational constraints we cannot study for 
large values of g (i.e., g > 2.25). 

We can set a cut-off for the anisotropy parameter A = 
Jiz/J\xy above which the essential physics for the system 
does not change much; for values of A above the cut-off, 
we expect that NN occupancy is in effect projected out. 
More precisely, we find that we can deal with large values 
of g (i.e., g > 2.25) by setting J\ z /J\ xv = Ao = 15 and 
still get the correct behavior of the observables thereby 
saving computational time. Fig.[4]shows m z versus h and 
S(n, 7 r) against m z for g = 2.5; results for different values 
of A are compared in order to decide on a cut-off value 
A = Ao. In the large anisotropic limit, a change in the 
anisotropy parameter dA can be shown to produce an 
additional effective field of h e g- = 2dA (at the mean-field 
level). So in Fig. [4](a), we shifted the h scale accordingly 
in an attempt to make the m z plots coincide. The good 
agreement between the two cases for A = 15 and A = 20 
gives us the freedom to use a cut-off of A = Ao = 15 at 
large values of g (i.e., g > 2.25). 

With the above simplification, we do the SSE-QMC 
simulation for larger values of g (such as g > 2.25) and 
observe phase-separated phases; we explain this based 
on Fig. [5] plotted at g = 2.5. At values of magnetic field 
h < 29.05, a single peaked structure occurs. On increas¬ 
ing h, at h = 29.062, a double-peaked structure results 
showing simultaneous existence of two phases (with mag¬ 
netizations centered at m z ~ .12 and m z ss .25). A fur¬ 
ther small increase to h = 29.07, leads to again a single 
peak (centered around m z ~ 0.26) signalling that a dis¬ 
continuous phase transition has occurred. Thus a phase 
separated state at g = 2.5 is clearly captured in Fig. [51(a). 


TABLE II. Autocorrelation times at g = 2.5 and t/iuo = 1. 


A = 10 

e = 1.0 

t = 2.5 

h = 18.0 

250289 

96613 

h = 18.5 

658170 

174995 

h = 19.0 

208295 

9835 

h = 19.5 

960 

113 


A = 15 

e = 4.0 

e = 6.0 

h = 28.0 

238406 

82608 

h = 28.5 

856138 

666353 

h = 29.0 

161336 

8646 

h = 29.1 

4847 

4778 


In Fig. [5](b), we show the evolution of the SF density ps 
and as h is varied and capture the first-order 

transition at h « 29.06. 

The results for the magnetization versus magnetic field 
in the transition region for values of g = 2.25, 2.5 and 3.0 
are shown in Figs. [61(b), (c), and (d). There is a gradual 
increase in the quantum of the magnetization jump as 
the coupling g is increased. We should also mention here 
that, particularly in the calculations involving large A, 
we wanted to ensure that the numerical data in adjacent 
bins are not correlated. To this end, we compute the 
autocorrelation time T int defined as2£ 

1 00 

T int\nriz\ = T ^ ' A m z jt), 
t=l 


where 


( t ) 


< m z (i + t)m z (i ) > — < m z (i) > 2 
< m z (i) 2 > — < m z (i) > 2 


(5) 


with i and t being Monte Carlo times defined in units of 
Monte Carlo sweeps (MCS). 

Typically, we see that for moderately high values of 
A, e = A/4 can restrict Tint from attaining very large 
values. Thus, for A < 10, we use a large bin size (i.e., 




FIG. 5. (Color online) First-order transition shown by 
the effective BH Hamiltonian when same-sublattice 
transport is inadequate, (a) Magnetization m z histograms 
in the transition region showing phase separation through a 
double-peaked structure at h = 29.062; (b) evolution of ps 
and S( 7r, 7r) vs h during phase transition. Both plots are for 
g = 2.5 and t/ujo = 1. 
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FIG. 6. (Color online) Evolution from continuous transi¬ 
tion to discontinuous transition as coupling strength 
increases in the BH model. Magnetization vs magnetic 
field in a 16 x 16 lattice for t/uio = 1 and for different values 
of g. 


15,00,000 MCS) in our simulations and keep Tint well 
within the bin size in order to produce meaningful re¬ 
sults from our simulation. However, for larger A, even 
choosing e = A/4 cannot keep autocorrelation times 
sufficiently smaller than such large bin sizes. So, in 
those cases, we use larger e values (i.e., e =6 and 8 for 
A=15 and 20 , respectively) and even larger bin sizes (i.e., 
22,00,000 MCS). The values of the autocorrelation times, 
for g = 2.5 and at various fields h close to the transition, 
are shown in Table El It should also be noted that we 
cannot take e too large, as a calculation with e larger 
than A does not produce meaningful results. 


Lastly, we mention that our phase diagrams, both 
for the XXZ model and our extension of it, are simi¬ 
lar (though not identical) to what were obtained using 
Lanczos method in Ref. [23]. 


V. SUMMARY 

In this work, we studied the effective Hamiltonian of a 
Bose-Holstein model using the SSE-QMC technique and 
obtained the ground state phase diagram. We found 
that supersolidity is realized at intermediate couplings; 
whereas, at large couplings the system phase separates 
because the double-hopping terms (that produce trans¬ 
port in the same sublattice) are not dominant. 

Our results on a large 16 x 16 lattice represent well the 
system behavior in the thermodynamic limit, as demon¬ 
strated in Fig. [T| using different finite-size calculations. 
Our results are only qualitatively similar to those ob¬ 
tained earlier using modified Lanczos technique on a 
much smaller 4x4 cluster—. 

We overcame computational difficulties for large repul¬ 
sive interactions by devising a cutoff repulsive strength; 
above the cutoff, the system properties (as shown in 
Fig. U) become essentially independent of the strength 
of repulsion (because repulsive interactions project out 
nearest-neighbor occupation). To mimic the results of 
statistically independent configurations, we considered 
sufficiently large e values and kept the auto-correlations 
within acceptable limits. 

Our work is an exercise in SSE-QMC study of a simple 
but important model which has significance in various 
fields. We hope that the present results will stimulate 
further investigations in allied areas such as frustrated 
quantum magnets in various geometries and at various 
magnetizations; coherence dynamics of excitons/spins in 
the presence of phonon environments (pertinent to quan¬ 
tum computation and artificial light harvesting); dimer 
formation and dimer correlations in Hubbard-Holstein 
models, HCBs coupled to multimode phonons^!, etc. 


Appendix A: SSE bond operators 

In SSE-QMC study of Heisenberg spin systems, the 
Hamiltonian is written as a bond Hamiltonian. Partic¬ 
ularly, in our case we write H = — X]i=i Sb where 
bi, 62 , & 63 denote the NN, NNN, and NNNN bonds in 
our spin model, respectively. Each of such Hi >i consists 
of the diagonal (H 1 ^) and the off-diagonal (L^.O parts 
and is given as = Hi ibi + with expressions 

H 1M =Ci- JizSfaSfa Y h b [Sf {bt) Y Sfa] 

H 2M =^[S+ bi) S- (bi) + H.c.], (Al) 

where J 2 Z — J 3 z — 0; C) = Jiz!^ Y ^b Y cJ\xyi e ^ 0, 
and hb = hJi xy /z with the coordination number z = 12. 
In our model, a two-spin matrix element of any of these 
operators can never become negative. 


1 S. H. Blanton, R. T. Collins, K. H. Kelleher, L. D. Rot- study of Bai x K x Bi03 from charge-density-wave insulator 
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Comparison between ED and SSE results 

We have benchmarked our QMC calculations by com¬ 
paring our calculated values with those obtained by 
exact-diagonalization (ED) methods. We consider both 
the XXZ model and its simple extension, namely, the 
J\ — J 2 model. We find that the energy, magnetization, 
structure factor S{ 7r,7r) and SF density ps of our SSE 
calculations match quite well with those from the ED re¬ 
sults. The comparisons of the calculated and ps 

for the Ji — J 2 model are shown in Fig. ISTT aVfbh 

Our SSE results also compare well with various world¬ 
line Monte Carlo results^— and also the SSE QMC 


results^ available in the literature. 
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FIG. SI. (Color online) Comparison of SSE with ED results 
on 4 x 4 clusters for a 2D Ji — J 2 model: (a) S(7r, 7 r) vs h (ED 
results are obtained using LAPACK), (b) ps vs J 2 /J 1 (ED 
results are taken from Refs, jl| and[§). 
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